| Variável | Descrição |
|---|---|
PropBaixaRenda |
Variável resposta binária (0 ou 1), derivada da “Proporção de pessoas com renda menor de ¼ do salário mínimo”. |
TxAnalfa |
Taxa de analfabetismo do município (%) |
TxDesemp |
Taxa de desemprego do município (%) |
PopRes |
População residente do município (Contagem) |
Saneat |
Quantidade de moradores com abastecimento de água (Contagem) |
ObtOcorr |
Óbitos ocorridos no município |
url <- "https://gist.githubusercontent.com/jonhsp/5cf24aa66c3fb173fea117bf4cb1dbc7/raw/29b6a4bcb197ab9be16663f8c73edcfd6b3c51f3/MunicipiosRS.xlsx%2520-%2520Quartis_1_e_3.csv"
df <- read_csv(url)
df %<>%
select(-'Prop<025SM') %>%
# Isolando o nome do municipio (remove codigo e espaco inicial)
mutate(MUNICIPIO = str_remove(MUNICIPIO, "^\\d+\\s")) %>%
mutate(MUNICIPIO = str_to_lower(MUNICIPIO)) %>%
# Tratando hifen como zero e convertendo para numerico
mutate(across(
.cols = c(TxAnalfa, TxDesemp, PopRes, ObtOcorr, Saneat),
.fns = ~ {
.x %>%
str_replace("^-$", "0") %>%
str_replace_all("\\.", "") %>%
str_replace_all(",", ".") %>%
as.numeric()
}
)) %>%
# Criando a variável binária como fator
mutate(y_bin = factor(PropBaixaRenda,
levels = c(0, 1),
labels = c("Baixa Proporção", "Alta Proporção"))) %>%
# Renomeando as colunas para minúsculo
rename_all(str_to_lower) %>%
# Selecionando variaveis
select(municipio, y_bin, txanalfa, txdesemp, popres, obtocorr, saneat)
summary(df)
## municipio y_bin txanalfa txdesemp
## Length:249 Baixa Proporção:124 Min. : 0.900 Min. : 0.0000
## Class :character Alta Proporção :125 1st Qu.: 3.200 1st Qu.: 0.8525
## Mode :character Median : 5.850 Median : 1.7150
## Mean : 6.645 Mean : 2.1223
## 3rd Qu.: 9.700 3rd Qu.: 2.9000
## Max. :19.100 Max. :10.9800
## NA's :1 NA's :1
## popres obtocorr saneat
## Min. : 0 Min. : 0.00 Min. : 0
## 1st Qu.: 2669 1st Qu.: 6.00 1st Qu.: 2656
## Median : 4171 Median : 13.00 Median : 4162
## Mean : 10238 Mean : 56.28 Mean : 10183
## 3rd Qu.: 8771 3rd Qu.: 49.00 3rd Qu.: 8747
## Max. :435564 Max. :2889.00 Max. :433266
##
df %>%
# Seleciona apenas colunas numéricas
select(where(is.numeric)) %>%
# Pivota os dados para o formato longo (ideal para facet_wrap)
pivot_longer(cols = everything(), names_to = "Variavel", values_to = "Valor") %>%
ggplot(aes(x = Valor)) +
geom_histogram(aes(y = ..density..), bins = 20, fill = "steelblue", color = "white", alpha = 0.8) +
geom_density(color = "red", size = 1) +
# Cria um mini-gráfico para cada variável
facet_wrap(~ Variavel, scales = "free") +
labs(
title = "Distribuição das Variáveis Preditoras",
x = "Valor da Variável",
y = "Densidade"
) +
theme_light()
df_transf <- df %>%
# Remove municípios com população zero (evita log(0) = -Inf)
filter(popres > 0) %>%
# Remoção de NAs
#drop_na() %>%
# log populacao
mutate(log_popres = log(popres)) %>%
# Obtos por 100.000 habitantes
mutate(obt_100k = obtocorr / popres * 100000) %>%
# Saneamento por 100.000 habitantes
mutate(saneat_100k = saneat / popres * 100000) %>%
# Remoção das variáveis originais
select(-c(popres, obtocorr, saneat))
df_transf %>%
# Seleciona apenas colunas numéricas
select(where(is.numeric)) %>%
# Pivota os dados para o formato longo (ideal para facet_wrap)
pivot_longer(cols = everything(), names_to = "Variavel", values_to = "Valor") %>%
ggplot(aes(x = Valor)) +
geom_histogram(aes(y = ..density..), bins = 20, fill = "steelblue", color = "white", alpha = 0.8) +
geom_density(color = "red", size = 1) +
# Cria um mini-gráfico para cada variável
facet_wrap(~ Variavel, scales = "free") +
labs(
title = "Distribuição das Variáveis Transformadas",
x = "Valor da Variável",
y = "Densidade"
) +
theme_light()
# Visualização das variáveis numéricas por nível da variável resposta
df_transf %>%
# 1. Seleciona numéricas E a variável resposta (y_bin)
select(y_bin, where(is.numeric)) %>%
# 2. Pivota apenas as numéricas, mantendo o y_bin fixo (cols = -y_bin)
pivot_longer(cols = -y_bin, names_to = "Variavel", values_to = "Valor") %>%
# 3. Adiciona 'fill = y_bin' para separar as cores por grupo
ggplot(aes(x = Valor, fill = y_bin, color = y_bin)) +
# position = "identity" permite que os histogramas se sobreponham (transparencia ajuda)
geom_histogram(aes(y = ..density..), bins = 20, alpha = 0.4, position = "identity") +
# Linha de densidade também separada por cor
geom_density(alpha = 0.2, size = 1) +
facet_wrap(~ Variavel, scales = "free") +
labs(
title = "Distribuição das Variáveis por Grupo (Baixa vs Alta Proporção)",
x = "Valor da Variável",
y = "Densidade",
fill = "Grupo",
color = "Grupo"
) +
theme_light() +
# Move a legenda para baixo para dar mais espaço aos gráficos
theme(legend.position = "bottom")
# Visualização das variáveis numéricas por nível da variável resposta (Boxplot)
df_transf %>%
# 1. Seleciona numéricas E a variável resposta (y_bin)
select(y_bin, where(is.numeric)) %>%
# 2. Pivota apenas as numéricas, mantendo o y_bin fixo
pivot_longer(cols = -y_bin, names_to = "Variavel", values_to = "Valor") %>%
# 3. MUDANÇA: Eixo X é o Grupo (y_bin), Eixo Y é o Valor numérico
ggplot(aes(x = y_bin, y = Valor, fill = y_bin)) +
# Geometria de Boxplot
geom_boxplot(alpha = 0.7, outlier.colour = "red", outlier.shape = 1) +
# Escalas livres no eixo Y para acomodar magnitudes diferentes
facet_wrap(~ Variavel, scales = "free_y") +
labs(
title = "Boxplots das Variáveis por Grupo",
x = "Grupo",
y = "Valor da Variável",
fill = "Grupo"
) +
theme_light() +
theme(legend.position = "bottom")
# Gráfico de Dispersão: Analfabetismo vs Desemprego por Grupo
df_transf %>%
ggplot(aes(x = txanalfa, y = txdesemp, color = y_bin)) +
# Adiciona os pontos (dispersão)
geom_point(alpha = 0.6, size = 2) +
labs(
title = "Relação entre Analfabetismo e Desemprego",
subtitle = "Segmentado pela Proporção de Baixa Renda (Y)",
x = "Taxa de Analfabetismo (%)",
y = "Taxa de Desemprego (%)",
color = "Situação (Y)"
) +
theme_light() +
theme(legend.position = "bottom")
Estamos construindo os modelos progressivamente. O objetivo é encontrar o equilíbrio entre simplicidade e capacidade explicativa.
O modelo nulo consiste na estimação apenas do \(\beta_{0}\) (intercepto), sem incluir os efeitos das variáveis preditoras. Este modelo será ajustado apenas para comparação da deviance residual.
ajuste_nulo <- glm(y_bin ~ 1, data = df_transf, family = binomial(link = "logit"))
summary(ajuste_nulo)
##
## Call:
## glm(formula = y_bin ~ 1, family = binomial(link = "logit"), data = df_transf)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.000 0.127 0 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 343.8 on 247 degrees of freedom
## Residual deviance: 343.8 on 247 degrees of freedom
## AIC: 345.8
##
## Number of Fisher Scoring iterations: 2
Este é o modelo base (saturado), assumindo que todas as relações são lineares.
# Ajuste do modelo apenas com efeitos principais lineares
ajuste_1 <- glm(y_bin ~ txanalfa + txdesemp + log_popres + saneat_100k + obt_100k,
data = df_transf,
family = binomial(link = "logit"))
# Sumário para ver coeficientes e significância individual
summary(ajuste_1)
##
## Call:
## glm(formula = y_bin ~ txanalfa + txdesemp + log_popres + saneat_100k +
## obt_100k, family = binomial(link = "logit"), data = df_transf)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 4.213e+01 1.019e+02 0.413 0.6792
## txanalfa 1.309e+00 2.008e-01 6.517 7.16e-11 ***
## txdesemp 7.262e-01 2.438e-01 2.979 0.0029 **
## log_popres 3.227e-01 4.258e-01 0.758 0.4485
## saneat_100k -5.364e-04 1.015e-03 -0.528 0.5973
## obt_100k -1.951e-03 1.231e-03 -1.585 0.1129
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 343.80 on 247 degrees of freedom
## Residual deviance: 89.62 on 242 degrees of freedom
## AIC: 101.62
##
## Number of Fisher Scoring iterations: 7
# Análise de Deviance (ANOVA) para ver a contribuição de cada variável ao ser adicionada
# O teste "Chisq" é o adequado para GLMs binomiais
anova(ajuste_1, test = "Chisq")
## Analysis of Deviance Table
##
## Model: binomial, link: logit
##
## Response: y_bin
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 247 343.80
## txanalfa 1 230.449 246 113.35 < 2.2e-16 ***
## txdesemp 1 20.541 245 92.81 5.838e-06 ***
## log_popres 1 0.037 244 92.77 0.84706
## saneat_100k 1 0.140 243 92.63 0.70779
## obt_100k 1 3.014 242 89.62 0.08257 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
A análise do modelo inicial (apenas efeitos principais lineares) revela os seguintes pontos:
txanalfa (Analfabetismo): Altamente
significativa (\(p < 0.001\)). O
coeficiente positivo (\(1.31\))
confirma que quanto maior o analfabetismo, maior a chance de o município
pertencer ao grupo de alta pobreza.txdesemp (Desemprego): Significativa
(\(p = 0.003\)). Coeficiente positivo
(\(0.73\)), indicando que o desemprego
também impulsiona a pobreza.log_popres (População): Não
significativa (\(p = 0.45\)). Isso
sugere que, em uma relação puramente linear, o tamanho do município não
distingue bem os grupos de renda. Isso reforça a necessidade de testar
termos quadráticos.saneat_100k e obt_100k:
Não apresentaram significância estatística neste modelo (\(p > 0.05\)).Conclusão: O modelo captura bem os efeitos de educação e trabalho, mas falha em capturar a influência da demografia (população) e infraestrutura (saneamento) em sua forma linear.
Com base na análise exploratória, adicionamos termos quadráticos para População e Desemprego, mantendo as demais variáveis lineares. O objetivo é testar se a relação com a pobreza é curvilínea.
# Ajuste do modelo com polinômios de 2º grau para as variáveis identificadas
ajuste_2 <- glm(y_bin ~ txanalfa + saneat_100k + obt_100k +
poly(txdesemp, 2) +
poly(log_popres, 2),
data = df_transf,
family = binomial(link = "logit"))
# Sumário para verificar a significância dos termos quadráticos
summary(ajuste_2)
##
## Call:
## glm(formula = y_bin ~ txanalfa + saneat_100k + obt_100k + poly(txdesemp,
## 2) + poly(log_popres, 2), family = binomial(link = "logit"),
## data = df_transf)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.809e+01 1.023e+02 0.275 0.78369
## txanalfa 1.332e+00 2.108e-01 6.321 2.6e-10 ***
## saneat_100k -3.436e-04 1.028e-03 -0.334 0.73817
## obt_100k -3.396e-03 1.667e-03 -2.037 0.04170 *
## poly(txdesemp, 2)1 4.882e+01 1.689e+01 2.891 0.00384 **
## poly(txdesemp, 2)2 3.683e+01 1.667e+01 2.209 0.02716 *
## poly(log_popres, 2)1 1.091e+00 8.316e+00 0.131 0.89559
## poly(log_popres, 2)2 -2.059e+01 9.309e+00 -2.211 0.02700 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 343.801 on 247 degrees of freedom
## Residual deviance: 80.534 on 240 degrees of freedom
## AIC: 96.534
##
## Number of Fisher Scoring iterations: 8
# Teste de Deviance (ANOVA) para comparar se o Modelo 2 é superior ao Modelo 1
# H0: Os modelos são iguais (a complexidade extra não ajuda)
# H1: O Modelo 2 é melhor
anova(ajuste_nulo, ajuste_1, ajuste_2, test = "Chisq")
## Analysis of Deviance Table
##
## Model 1: y_bin ~ 1
## Model 2: y_bin ~ txanalfa + txdesemp + log_popres + saneat_100k + obt_100k
## Model 3: y_bin ~ txanalfa + saneat_100k + obt_100k + poly(txdesemp, 2) +
## poly(log_popres, 2)
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 247 343.80
## 2 242 89.62 5 254.181 < 2e-16 ***
## 3 240 80.53 2 9.086 0.01064 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Os resultados obtidos a partir do ajuste do modelo com termos quadráticos indicam uma superioridade clara sobre o modelo linear simples (Modelo 1). A análise detalhada segue abaixo:
Qualidade do Ajuste:
Teste de Razão de Verossimilhança (ANOVA): A comparação formal entre os modelos resultou em um p-valor de 0.0106 (Chi-quadrado). Como \(p < 0.05\), rejeitamos a hipótese nula de que os modelos são equivalentes. Isso prova estatisticamente que a inclusão da curvatura para Desemprego e População adicionou poder explicativo real ao modelo, e não apenas ruído.
Análise dos Coeficientes:
poly(txdesemp, 2)2, p=0.027) e População
(poly(log_popres, 2)2, p=0.027) foram estatisticamente
significativos. Isso sugere que a relação dessas variáveis com a pobreza
não é constante, variando conforme a magnitude dos indicadores (formato
de curva).txanalfa) permanece como o preditor mais robusto do modelo
(\(p < 0.001\)), mantendo seu efeito
positivo forte.saneat_100k) apresentou p-valor alto (0.73), sugerindo
que, na presença das outras variáveis (especialmente Analfabetismo), ela
não adiciona informação nova relevante para a classificação.obt_100k) tornou-se significativa (\(p=0.04\)) neste modelo ajustado.Conclusão: O Modelo 2 será adotado como base para a seleção final. O próximo passo será refinar este modelo aplicando o método Stepwise para remover automaticamente os preditores não significativos (como o Saneamento) e obter o modelo final parcimonioso.
Nesta etapa, refinamos o modelo quadrático removendo variáveis redundantes pelo critério de Akaike (AIC) e realizamos os diagnósticos de qualidade do ajuste.
Utilizamos o método stepwise (direção “both”) para remover preditores não significativos e chegar ao modelo mais parcimonioso possível.
# Aplica o algoritmo Stepwise no Modelo 2 (Quadrático)
# O R testará remover/adicionar variáveis para baixar o AIC
modelo_step <- step(ajuste_2, direction = "both", trace = 1)
## Start: AIC=96.53
## y_bin ~ txanalfa + saneat_100k + obt_100k + poly(txdesemp, 2) +
## poly(log_popres, 2)
## Df Deviance AIC
## - saneat_100k 1 80.645 94.645
## <none> 80.534 96.534
## - poly(log_popres, 2) 2 86.426 98.426
## - obt_100k 1 86.482 100.482
## - poly(txdesemp, 2) 2 101.088 113.088
## - txanalfa 1 276.631 290.631
##
## Step: AIC=94.65
## y_bin ~ txanalfa + obt_100k + poly(txdesemp, 2) + poly(log_popres,
## 2)
## Df Deviance AIC
## <none> 80.645 94.645
## + saneat_100k 1 80.534 96.534
## - poly(log_popres, 2) 2 86.943 96.943
## - obt_100k 1 86.579 98.579
## - poly(txdesemp, 2) 2 101.596 111.596
## - txanalfa 1 277.305 289.305
# Exibe o resumo do modelo vencedor
summary(modelo_step)
##
## Call:
## glm(formula = y_bin ~ txanalfa + obt_100k + poly(txdesemp, 2) +
## poly(log_popres, 2), family = binomial(link = "logit"), data = df_transf)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -6.124315 1.168673 -5.240 1.60e-07 ***
## txanalfa 1.334196 0.210857 6.327 2.49e-10 ***
## obt_100k -0.003382 0.001666 -2.030 0.04237 *
## poly(txdesemp, 2)1 49.542674 16.900795 2.931 0.00337 **
## poly(txdesemp, 2)2 37.308868 16.737594 2.229 0.02581 *
## poly(log_popres, 2)1 1.500561 8.255232 0.182 0.85576
## poly(log_popres, 2)2 -20.826304 9.307757 -2.238 0.02525 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 343.801 on 247 degrees of freedom
## Residual deviance: 80.645 on 241 degrees of freedom
## AIC: 94.645
##
## Number of Fisher Scoring iterations: 8
O modelo capta a não linearidade dignificativa de log populção e da taxa de desemprego, porém também acrescenta a taxa de obtos ao nível de 5% de significancia, porém a mesma não seria incluida ao nível de 1% de significancia. O \(\beta_{obtos}\) foi estimado em -0.003, o que indica que a taxa de obtos tem um efeito negativo na probabilidade de alta taxa proporção de pobreza. Esse resultado calsa estranheza, pois indica que municípios mais violetos tendem a ter menores taxas de pobreza, contudo, a significância desta variável deve estar relacionada com o efeito de outros fatores, externos ao modelo, que parcialmente são representandos pela quantidade de obtos.
Visando a simplicidade do modelo e a interpretabilidade dos
resultados, a taxa de obtos será removida, sem grandes perdas ao modelo,
visto que (AIC) < 4
modelo_final <- update(modelo_step, . ~ . - obt_100k)
# Exibe o resumo do modelo vencedor
cat("----Summary do Modelo Final----")
## ----Summary do Modelo Final----
summary(modelo_final)
##
## Call:
## glm(formula = y_bin ~ txanalfa + poly(txdesemp, 2) + poly(log_popres,
## 2), family = binomial(link = "logit"), data = df_transf)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -7.1293 1.1086 -6.431 1.27e-10 ***
## txanalfa 1.2417 0.1895 6.553 5.66e-11 ***
## poly(txdesemp, 2)1 42.5755 15.1374 2.813 0.00491 **
## poly(txdesemp, 2)2 29.5031 14.6302 2.017 0.04374 *
## poly(log_popres, 2)1 -5.4882 8.1155 -0.676 0.49888
## poly(log_popres, 2)2 -13.3643 8.9366 -1.495 0.13479
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 343.801 on 247 degrees of freedom
## Residual deviance: 86.579 on 242 degrees of freedom
## AIC: 98.579
##
## Number of Fisher Scoring iterations: 8
cat("----Anova entre os Modelos----")
## ----Anova entre os Modelos----
anova(modelo_final,modelo_step, ajuste_2, test = "Chisq")
## Analysis of Deviance Table
##
## Model 1: y_bin ~ txanalfa + poly(txdesemp, 2) + poly(log_popres, 2)
## Model 2: y_bin ~ txanalfa + obt_100k + poly(txdesemp, 2) + poly(log_popres,
## 2)
## Model 3: y_bin ~ txanalfa + saneat_100k + obt_100k + poly(txdesemp, 2) +
## poly(log_popres, 2)
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 242 86.579
## 2 241 80.645 1 5.9332 0.01486 *
## 3 240 80.534 1 0.1115 0.73850
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
cat("---- AIC ----")
## ---- AIC ----
AIC(ajuste_nulo, ajuste_1, ajuste_2, modelo_step, modelo_final)
## df AIC
## ajuste_nulo 1 345.80100
## ajuste_1 6 101.62016
## ajuste_2 8 96.53401
## modelo_step 7 94.64546
## modelo_final 6 98.57871
A tabela abaixo apresenta a Razão de Chances (Odds Ratios) para as variáveis do modelo final.
# Extração dos coeficientes e cálculo do OR (exp(beta))
tabela_or <- exp(cbind(OR = coef(modelo_final), confint(modelo_final)))
# Exibição formatada
knitr::kable(tabela_or, digits = 3, caption = "Razão de Chances (Odds Ratios) e IC 95%")
| OR | 2.5 % | 97.5 % | |
|---|---|---|---|
| (Intercept) | 1.000000e-03 | 0.000 | 5.000000e-03 |
| txanalfa | 3.462000e+00 | 2.508 | 5.334000e+00 |
| poly(txdesemp, 2)1 | 3.092593e+18 | 55970162.875 | 1.107302e+33 |
| poly(txdesemp, 2)2 | 6.501656e+12 | 143.121 | 4.430122e+26 |
| poly(log_popres, 2)1 | 4.000000e-03 | 0.000 | 2.614928e+04 |
| poly(log_popres, 2)2 | 0.000000e+00 | 0.000 | 1.210000e+01 |
# 1. Efeito do Analfabetismo (Linear)
# Mostra a probabilidade predita de Alta Pobreza conforme a Taxa de Analfabetismo aumenta
plot_analfa <- ggpredict(modelo_final, terms = "txanalfa") %>%
plot() +
labs(
title = "Efeito do Analfabetismo",
x = "Taxa de Analfabetismo (%)",
y = "Probabilidade de Alta Pobreza"
) +
theme_light()
# 2. Efeito do Desemprego (Quadrático)
# Mostra a curva não-linear capturada pelo polinômio
plot_desemp <- ggpredict(modelo_final, terms = "txdesemp [all]") %>%
plot() +
labs(
title = "Efeito do Desemprego (Não-linear)",
x = "Taxa de Desemprego (%)",
y = "Probabilidade de Alta Pobreza"
) +
theme_light()
# 3. Efeito da População (Quadrático)
# Plota a curvatura ajustada para o Log da População
plot_pop <- ggpredict(modelo_final, terms = "log_popres [all]") %>%
plot() +
labs(
title = "Efeito do Tamanho da População",
x = "Log(População)",
y = "Probabilidade de Alta Pobreza"
) +
theme_light()
# Exibir os gráficos juntos usando gridExtra
# Organiza em 2 colunas: Analfabetismo e Desemprego na linha 1, População centralizada na linha 2
grid.arrange(plot_analfa, plot_desemp, plot_pop, ncol = 3,
layout_matrix = rbind(c(1, 2), c(3, 3)))
# 1. Efeito do Analfabetismo (Linear)
# Mostra a probabilidade predita de Alta Pobreza conforme a Taxa de Analfabetismo aumenta
plot_analfa <- ggpredict(modelo_final, terms = "txanalfa", typical = 'mode') %>%
plot() +
labs(
title = "Efeito do Analfabetismo",
x = "Taxa de Analfabetismo (%)",
y = "Probabilidade de Alta Pobreza"
) +
theme_light()
# 2. Efeito do Desemprego (Quadrático)
# Mostra a curva não-linear capturada pelo polinômio
plot_desemp <- ggpredict(modelo_final, terms = "txdesemp [all]", typical = 'mode') %>%
plot() +
labs(
title = "Efeito do Desemprego (Não-linear)",
x = "Taxa de Desemprego (%)",
y = "Probabilidade de Alta Pobreza"
) +
theme_light()
# 3. Efeito da População (Quadrático)
# Plota a curvatura ajustada para o Log da População
plot_pop <- ggpredict(modelo_final, terms = "log_popres [all]", typical = 'mode') %>%
plot() +
labs(
title = "Efeito do Tamanho da População",
x = "Log(População)",
y = "Probabilidade de Alta Pobreza"
) +
theme_light()
# Exibir os gráficos juntos usando gridExtra
# Organiza em 2 colunas: Analfabetismo e Desemprego na linha 1, População centralizada na linha 2
grid.arrange(plot_analfa, plot_desemp, plot_pop, ncol = 3,
layout_matrix = rbind(c(1, 2), c(3, 3)))
# 1. Preparar os dados com as predições
# Usamos os dados reais para ver onde os municípios se posicionam no espaço 3D
df_3d <- df_transf %>%
mutate(
prob_predita = predict(modelo_final, type = "response"),
# Criamos um texto formatado para quando passar o mouse
hover_text = paste(
"<b>Município:</b>", municipio,
"<br><b>Analfabetismo:</b>", txanalfa, "%",
"<br><b>Desemprego:</b>", txdesemp, "%",
"<br><b>População:</b>", exp(log_popres) %>% round(0), # Revertendo o log para ler fácil
"<br><b>Prob. Alta Pobreza:</b>", scales::percent(prob_predita, accuracy = 0.1)
)
)
# 2. Gerar o Gráfico de Dispersão 3D
plot_ly(data = df_3d,
x = ~txanalfa,
y = ~txdesemp,
z = ~log_popres,
color = ~prob_predita,
colors = c("#440154", "#21908C", "#FDE725"), # Escala Viridis (escura -> clara)
type = "scatter3d",
mode = "markers",
marker = list(
size = 6, # Tamanho das bolinhas
opacity = 0.8, # Transparência para ver pontos atrás
line = list(width = 0) # Remove borda para ficar mais limpo
),
text = ~hover_text,
hoverinfo = "text"
) %>%
layout(
title = "Espaço Tridimensional dos Preditores<br><sub>Cor = Probabilidade de Alta Pobreza</sub>",
scene = list(
xaxis = list(title = "Analfabetismo (%)"),
yaxis = list(title = "Desemprego (%)"),
zaxis = list(title = "Log(População)"),
camera = list(eye = list(x = 1.5, y = 1.5, z = 1.2)) # Posição inicial da câmera
)
)
# 1. Preparar os dados de classificação
df_classificacao <- df_transf %>%
mutate(
# Probabilidade predita pelo modelo
prob = predict(modelo_final, type = "response"),
# Classificação baseada no limiar de 50% (0.5)
# Se prob > 0.5 -> Classifica como "Alta Proporção" (nível 2 do fator)
# Se prob <= 0.5 -> Classifica como "Baixa Proporção" (nível 1 do fator)
classe_predita = ifelse(prob > 0.5, "Alta Proporção", "Baixa Proporção"),
# Verificar se acertou ou errou
# Como y_bin é fator, convertemos para caracter para comparar
resultado = ifelse(as.character(y_bin) == classe_predita, "Acerto", "Erro"),
# Detalhe do erro (Falso Positivo vs Falso Negativo) para o tooltip
tipo_erro = case_when(
resultado == "Acerto" ~ "Classificação Correta",
classe_predita == "Alta Proporção" ~ "Falso Positivo (Predisse Alta, mas é Baixa)",
TRUE ~ "Falso Negativo (Predisse Baixa, mas é Alta)"
),
# Texto para o mouse
hover_text = paste(
"<b>Município:</b>", municipio,
"<br><b>Analfabetismo:</b>", txanalfa, "%",
"<br><b>Desemprego:</b>", txdesemp, "%",
"<br><b>População:</b>", exp(log_popres) %>% round(0), # Revertendo o log para ler fácil
"<br><b>Prob. Alta Pobreza:</b>", scales::percent(prob, accuracy = 0.1)
)
)
# 2. Gerar o Gráfico 3D
plot_ly(data = df_classificacao,
x = ~txanalfa,
y = ~txdesemp,
z = ~log_popres,
color = ~resultado,
colors = c("#1f77b4", "#d62728"), # Azul para Acerto, Vermelho para Erro
symbol = ~resultado,
symbols = c("circle", "circle"), # Bolinha para Acerto, X para Erro
type = "scatter3d",
mode = "markers",
marker = list(
size = 6,
opacity = 0.8,
line = list(width = 1, color = "black") # Borda para destacar
),
text = ~hover_text,
hoverinfo = "text"
) %>%
layout(
title = "Diagnóstico de Classificação 3D<br><sub>Pontos vermelhos indicam onde o modelo errou</sub>",
scene = list(
xaxis = list(title = "Analfabetismo (%)"),
yaxis = list(title = "Desemprego (%)"),
zaxis = list(title = "Log(População)")
),
legend = list(title = list(text = "Resultado"))
)
# Definindo as medianas para fixar as variáveis de controle
med_analfa <- median(df_transf$txanalfa)
med_desemp <- median(df_transf$txdesemp)
med_pop <- median(df_transf$log_popres)
# ----------------------------------------------------------------------------
# GRÁFICO 1: Analfabetismo vs Desemprego (Fixando População)
# ----------------------------------------------------------------------------
grid1 <- expand.grid(
txanalfa = seq(min(df_transf$txanalfa), max(df_transf$txanalfa), length.out = 50),
txdesemp = seq(min(df_transf$txdesemp), max(df_transf$txdesemp), length.out = 50),
log_popres = med_pop
)
grid1$prob <- predict(modelo_final, newdata = grid1, type = "response")
p1 <- ggplot(grid1, aes(x = txanalfa, y = txdesemp)) +
geom_tile(aes(fill = prob)) +
geom_contour(aes(z = prob), color = "white", breaks = 0.5, size=0.8) +
geom_point(data = df_transf, aes(color = y_bin), alpha = 0.3, size = 1.5) + # Pontos reais
scale_fill_gradient2(low = "#2c7bb6", mid = "#ffffbf", high = "#d7191c", midpoint = 0.5, limits=c(0,1)) +
scale_color_manual(values = c("black", "darkred")) +
labs(title = "Analfabetismo vs Desemprego", subtitle = "População fixa na mediana") +
theme_minimal() + theme(legend.position = "none")
# ----------------------------------------------------------------------------
# GRÁFICO 2: Analfabetismo vs Log(População) (Fixando Desemprego)
# ----------------------------------------------------------------------------
grid2 <- expand.grid(
txanalfa = seq(min(df_transf$txanalfa), max(df_transf$txanalfa), length.out = 50),
log_popres = seq(min(df_transf$log_popres), max(df_transf$log_popres), length.out = 50),
txdesemp = med_desemp
)
grid2$prob <- predict(modelo_final, newdata = grid2, type = "response")
p2 <- ggplot(grid2, aes(x = txanalfa, y = log_popres)) +
geom_tile(aes(fill = prob)) +
geom_contour(aes(z = prob), color = "white", breaks = 0.5, size=0.8) +
geom_point(data = df_transf, aes(color = y_bin), alpha = 0.3, size = 1.5) +
scale_fill_gradient2(low = "#2c7bb6", mid = "#ffffbf", high = "#d7191c", midpoint = 0.5, limits=c(0,1)) +
scale_color_manual(values = c("black", "darkred")) +
labs(title = "Analfabetismo vs População", subtitle = "Desemprego fixo na mediana") +
theme_minimal() + theme(legend.position = "none")
# ----------------------------------------------------------------------------
# GRÁFICO 3: Desemprego vs Log(População) (Fixando Analfabetismo)
# ----------------------------------------------------------------------------
grid3 <- expand.grid(
txdesemp = seq(min(df_transf$txdesemp), max(df_transf$txdesemp), length.out = 50),
log_popres = seq(min(df_transf$log_popres), max(df_transf$log_popres), length.out = 50),
txanalfa = med_analfa
)
grid3$prob <- predict(modelo_final, newdata = grid3, type = "response")
p3 <- ggplot(grid3, aes(x = txdesemp, y = log_popres)) +
geom_tile(aes(fill = prob)) +
geom_contour(aes(z = prob), color = "white", breaks = 0.5, size=0.8) +
geom_point(data = df_transf, aes(color = y_bin), alpha = 0.3, size = 1.5) +
scale_fill_gradient2(low = "#2c7bb6", mid = "#ffffbf", high = "#d7191c", midpoint = 0.5, limits=c(0,1), name="Prob.\nPredita") +
scale_color_manual(values = c("black", "darkred"), name="Real") +
labs(title = "Desemprego vs População", subtitle = "Analfabetismo fixo na mediana") +
theme_minimal() + theme(legend.position = "right")
# ----------------------------------------------------------------------------
# Exibição Conjunta
# ----------------------------------------------------------------------------
# Organiza os 3 plots. O p3 fornece a legenda para todos.
grid.arrange(p1, p2, p3, ncol = 3)
car::vif(modelo_final)
## GVIF Df GVIF^(1/(2*Df))
## txanalfa 1.217938 1 1.103602
## poly(txdesemp, 2) 1.697873 2 1.141501
## poly(log_popres, 2) 1.591663 2 1.123215
shapiro.test(residuals(modelo_final))
##
## Shapiro-Wilk normality test
##
## data: residuals(modelo_final)
## W = 0.84722, p-value = 6.342e-15
O GVIF, bem como o teste de Shapiro Wilk, confirmam que não há evidência de multicolinearidade entre as variáveis preditoras.
hnp(modelo_final, resid.type = "deviance", how.many.out = TRUE, paint.out = TRUE,
main = "Envelope Simulado", xlab = "Percentil Teórico", ylab = "Resíduos Deviance")
## Binomial model
## Total points: 248
## Points out of envelope: 0 ( 0 %)
plot(residuals(modelo_final))